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1. Introduction and motivation 



A realistic setup for studying many non-perturbative QCD properties is obtained by introduc- 
ing two quark pairs, a light (/) mass degenerate one (u and d flavours: 2) and a heavier (/?) mass 
non-degenerate one (s and c flavours: 1 + 1): for short Nf = 2 + 1 + 1. Monte Carlo simulations in 
this setup are becoming increasingly feasible in the tmLQCD formulation [[[J ||, ||] with action 

5 = S G [U) + % Dj[U] Yi + Yh D h [U] Yh , V/ = [M] , Wh = M . (1.1) 
Here Sq is a suitable pure gauge action , the Dirac operators read (in the "physical" quark basis) 



D, = yV - iy$TiW cr + Hi , D h = yV - iy 5 T\ W cr + jj. h - £ A t 3 



W cr 



--V*V+M cr , 



(1.2) 
(1.3) 



and £/, is the bare s-c quark mass splitting. Among the general features of the above formulation [^] 
we recall: automatic 0(a) improvement robust quark mass protection against "exceptional con- 
figurations" [jj]], expected moderate CPU-cost for unquenched simulations (assuming metastability 
problems related to the lattice phase structure are solved) IgJ. The determinant of Df, (eq. (|L~2|)) is 
real and positive provided |e/,| < |jU/,|, which, if Z P /Z S < 1, induces some limitation on the renor- 
malised quark masses, m C}S = Zp 1 ^ ± ^£/j), one can simulate. 

Through rescaling and chiral rotations of the quark fields (which do not affect the fermion 
determinant, besides an irrelevant constant factor) and by setting fli = 2k„\Ii ./, , 6/^ = 2K cr Eij 1 , the 
lattice Dirac operators (eq. (1.2)) can be rewritten in a form more convenient for MC simulations: 

(Di) 2x 2 = [yV + W cr ] 2K- r + //i/7 5 T3 , {D h ) 2x2 = [yV + W CT ] 2fc cr + ^75^ + £/,Ti . (1.4) 

As the value of K cr is a priori unknown, K cr is generically replaced by K in the definition of fli and 
£/,/! 2 - With these premises, we focus on the 75-Hermitian partner of the Dirac operator(s) above, 



22x2 







Q+ £75 


£7s Q ~ ifi 




£7s Q 



7 V_^V*V + ^-4 



where we have for the moment dropped the quark pair labels / and h. In the mass degenerate case 
(£ = 0), a standard HMC algorithm can be employed in view of 

det[e 2 x 2 ] = det[<2+<2-] = det[(& + p 2 ] ^ J ^ e -^{Q + Q-TH j 

with (p a single flavour pseudofermion field. This fits well our needs for the I quark pair. In the 
mass non-degenerate case (£ / 0) no plain HMC algorithm is straightforwardly applicable because 

det[(2 2 x 2 ] =det[<2+75<2--£ 2 75] = det[e+2- - ^G+T^Vs] , (1.6) 

cannot be reproduced by / with A a 1 -flavour matrix. Owing to the flavour non- 

diagonal structure of the matrix (L5), we propose to deal with a two-flavour matrix, 22x2^2x2' 
acting on a two-flavour pseudofermion field, <I>, and to use a polynomial P(s) ~ 1/ sfs, which gives 



det[(2 2 x2] O / ^<J> e-*^^^- 



4> 



(1.7) 



1 In unquenched computations the choice of S g [U ] is important for the phase structure of the lattice model [H, 
2 Then tmLQCD at maximal twist is obtained by setting K to a sensible estimate of Kc r for all quark pairs pp. 
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at least in the <I>-heatbath and the reweighting (or acceptance) correction. A PHMC algorithm 
appears thus a natural choice to include the effects of the h quark pair in MC simulations. 

2. A preconditioned (P)HMC algorithm for N f = 2 + 1 + 1 flavours 

We outline here the mixed HMC-PHMC (briefly (P)HMC) algorithm we are developing, which 
incorporates even-odd (EO) and mass-shift preconditionings Using different Molecular Dy- 
namics (MD) time steps for different MD force contributions we expect to obtain a good perfor- 
mance, in line with that recently achieved in simulations with two Wilson quark flavours || \R 

Denoting by IT the momenta conjugated to the gauge field, the MD-Hamiltonian reads 

1. 



H 



2+1 + 1 



-n n+s G +^lP [Q h Q h ] v h 



where P(Qi,Qt,) is a polynomial in QhQ\ approximating [<2/i<2jj ^ 2 (see sect. |J) while in the I quark 



sector mass-shift preconditioning is applied \q, |10J on top of EO preconditioning (see eq. ( |2.4[ ) and 
eq. ( |2.5[ )). In this way we need two pseudofermions, </>i and 02 > for the / quark sector and a two- 
flavour pseudofermion field, <!>/,, for the h quark sector. The EO preconditioned Dirac operators are 



Qh = j5 



h ( i + 



MoeM eo 

M oc (l+i/j;,y 5 )M et 



(M eo(oe) )^ = -k£ [(i + y M )^ Cy)^-/i + (i - y^u^dy^p] , 



with 2/, having a 2 x 2 flavour structure that is made apparent in the r.h.s. of eq. ( |2.2[ ), and 

M oe (l-i(fii + 8fii)Y 5 )M e , 



Q'l = 7s 



l + i(ft + 5ft)%- 



Q'l = 75 



l+i'A/75 



l + (ft + 5ft) 2 
M oe (\-ijxiy 5 )M eo 



(2.2) 
(2.3) 

(2.4) 
(2.5) 



with £2/ and £2" carrying the shifted (£/ + djli) and the physical (£/) twisted mass parameters. We 
remark that (due to the absence of the Sheikholeslami-Wohlert term in the fermionic action) the 



gauge field enters the Dirac matrices eqs. ( |2.2| ), ( |2.4| ) and ( |2.5| ) only through M eo and M oe , eq. (|2.3|). 
It follows that the evaluation of the MD driving force fl = — 5f/#2+i+i can be, as usual, traced back 
to that of 8uM eo and 8(jM oe plus the (many) necessary applications of the relevant Dirac matrices. 



With the (P)HMC update [0, [10(] dictated by the Hamiltonian #2+1+1 (eq. fl2.1|)) one ends (af- 
ter thermalisation) with a sample of gauge configurations equilibrated with respect to the effective 
gauge action Sc[U] — log | det<2'/| 2 + logdetP(<2/i<2/[)- A way to correct for the polynomial approx- 
imation of [QhQ' h ]~^ 2 is to reweight all the observables G = ff\U] with a correction factor [[J] that 
provides a noisy estimate of det[QhQl]^ 2 P(QhQl) = &&QhP{QhQ\) (see also sect. |3|). 
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2.1 MD force contributions and multiple time scales 

The MD driving force can be denned (omitting for brevity all indices) as 

n = -8uH 2 + w =tr [5I/F + FW] , F = F G + F h + F n +F a , (2.6) 



where (see eq. (2.1)) the pure gauge contribution {Fq) is completely standard, the / quark sector 



contributions (Fn and Fa) can be straightforwardly evaluated following Refs. [|10j, |11J] and also the 
h quark sector contribution (F/,) poses no principle problems (see sect. || for more details). 

Based on Refs. []9| [HS, we expect that for typical choices of So, the lattice spacing and the 
quark masses, one can have a hierarchy in the average (av) size of the individual force contributions, 

\F G \av > \Fll\av > \Fn\av , \Fl\av > \Fh\av > \Fn\av , (2-7) 

provided fa > is not too small and 8 fa > is appropriately chosen [JTot] - \Fh\ av is expected to be 
small for heavy quark masses, i.e. large \pLh\ and \fa t \ > |8/,| 3 . 



The hierarchy in eq. ( |2.7| ) -once realised- suggests that an optimal performance is obtained by 
implementing a MD leapfrog scheme where the force contributions (Fg,F/,,F/i,F/2) enter associ- 
ated to different time steps (Stq, 8%, Stu, 8112) so as to get 

\F G \av8T G - \F h \av8Th ^ \Fn\av8z n ^ \Fa\ m Sx a , SxiNi = T traj , i = {G,ll,l2,h} , (2.8) 

with (Nc,Nh,Nn,Ni2) a set of integers and T tra j ~ 1 the time length of a MD trajectory. 

2.2 Some possible algorithmic variants 

Several modifications of the above (P)HMC algorithmic scheme are of course possible. For 
instance, the correction for the polynomial approximation can be moved, fully or partially, from 
the reweighting to a modified A/R Metropolis step Jl2| , |l3| ]. If this is done "fully" the modified 



A/R step compensates for both det[y QhQlP(QhQl)] / 1 and the finite MD time step(s). 

The update of the gauge field itself can be performed by means of a Multiboson-like algo- 



rithm [14] using suitable polynomials to approximate the appropriate power of the inverse Her- 



mitean Dirac matrices (see e.g. Ref. [13] for more details) relevant for the / and h quark sectors. 

Another interesting possibility is to employ a non-standard HMC algorithm, whose MD is 
guided by an Hamiltonian H 7^ Z/2+1+1 sucn that a good acceptance is still obtained in the A/R 



test. One might try e.g. an H that differs from flj+i+i (eq. (|2JJ)) only by the replacement of 
&l P(QhQh)h^h with Qj^^u- Numerical experience is of course crucial to test these possible 
variants and choose the most efficient one. 



3. Polynomial approximations for the h quark sector (s and c flavours) 

The well known Chebyshev polynomial approximation method allows to approximate the in- 
verse square root of the operator S = PhQhQt> where p/, is a positive normalisation factor such that 
on "practically all" gauge configurations the highest eigenvalue of S, say sh, satisfies 0.8 < Sh < 1. 

3 It is important that |e;,| is not too close to \ph\. In fact, for |e;,| > not only the positivity of the determinant is 
no longer guaranteed, but the matrix Qi, (as well as £)/,) can even develop zero eigenvalues. 
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The polynomial P(5) = P njSL (S) of even degree n in 5 that is designed to approximate S 1 in the 
eigenvalue interval [sl, 1] 4 can be written (via a product representation with a normalisation factor 
,jY > and roots zu, k = l,...,n such that z n +i-k = z* k ) in a manifestly positive form 



/t=i 



B n/2 >SL B n /2, SL 



S-^ 2 [l+R, hSL ] , (3.1) 



where /?„ iL = R nSL (S) is called the relative fit error. Denoting by s a generic eigenvalue of 5, 
in Fig. H (left panel) we plot for illustration the relative fit error R IUSL (s) for si = 10~ 4 and a 
value of n, namely n = 162, taken such that \R„ , SL (s)\ < 0.03. The chosen value of Sl is rather 
conservative, since it is very low as compared to those that, based on simulation experience with 
two mass degenerate quarks, we expect to have to face in the h quark sector while working with 
realistic parameters. We see from Fig. [j] that R, h!iL (s) tends to increase when decreasing s £ [sl, 1], 
which we believe is acceptable in view of the expected non-high density of eigenvalues in the 
low end of the spectrum of S. A conservative and an "effective" measure of the magnitude of 
the relative fit error are thus given, respectively, by 8m = vaax se t SL u \R, USL (s)\ = \R n ,s L ( s L)\ an d 
Suv = max .«e[o.5,i] l^n,s L (' ? )l with Suv smaller than Sir (typically by an order of magnitude). 



0.005 



■0.01 



n = 162 






s L =10 4 

imwmimmmtiitm 




5, R = 3.0 10" 2 
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0.4 0.6 
s 



or 
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0.001 0.002 0.003 0.004 
S 



0.005 



Figure 1: R(s) = R n . SL (s) for (n,s L ) = (162, 1(T 4 ) (left) and R(s) = R H ,s L (s) for (n,s L ) = (620, 10" 4 ) (right). 

For the task of computing the contribution Fy, to the MD driving force we plan to exploit the 
product representation of P n SL (S), see eq. (|Q|), and then proceed analogously to Ref. A careful 



ordering [15] of the n monomials in S together with 64 bit precision should be sufficient to keep 
rounding errors under control for polynomials with n up to one thousand. 



There are two places where a second polynomial approximation (or an equivalent method such 



as a rational CG solver [ 16]) is needed. The first place is the generation of <!>/, distributed according 
to exp[— &jP ntSL (§)<&f l ]> see ec l- If 77, is a random Gaussian (two flavour) vector, we can 

construct <!>/, = B~L , (5) n, by e.g. evaluating 



n/2,si 
it 



-1 



> h = P, lSL (S)Bl /2sL (S) \^p- h Q h ] r h , P n , L (S) = [VSP n . SL (S)\ [l +R~n, L (S)] , (3.2) 

where Phj l (§) is a new polynomial of degree n in S that approximates [VsP n ^ L (S)]~ l with very 
high precision in the eigenvalue interval [sl, 1]. Values of the corresponding relative fit error, R~n )SL , 

4 By si we denote the lowest eigenvalue of S. It can be estimated in the early stages of a simulation, e.g. by starting 
with a "trial" polynomial P N ^ and then changing N and A on the basis of on-the-fly measurements of eigenvalues of S. 
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not exceeding in modulus 5jr = 10~ 6 can be reached for a degree n = 620, in the above mentioned 
case of sl = 10~ 4 and n = 162, as it is illustrated in Fig. [I] (right panel). The evaluation of P n ^ L {S) 
times a vector can be safely performed in 64 bit arithmetics by using a recursive relation (such as 
Clenshaw's) that is stable against roundoff. 

The other place where the high precision polynomial P thSL (S) in eq. (3.2) may be useful is in 
the evaluation of W [j\ ; 17] = exp j 7] t (l - [P„ iSL (5) Vs] ~ 1 ) rj } , the noisy reweighting factor that is 
needed to get the "exact" v.e.v. of a generic observable G = &[U] (here Y\ is the noise field) [Q] 

(0)=fDn[U]fDri e-^n W[ri;U] 0]U] / fDn[U]fDr} W[r);U] , 

jD^[U]=jDUe- s ^ d e t[Q';Qf[U]]d e t[P(S[U])]' 1 = Z {P)HUC . (3.3) 

4. Conclusions and Acknowledgements 

We discussed an exact algorithm for Nf ■ = 2 + 1 + 1 flavours of maximally twisted quarks. 
Implementation of the algorithm and investigation of important numerical properties, such as the 
spectrum of QhQ\ an d the magnitude of the contribution Ff, to the MD driving force, are in progress. 

We thank K. Jansen for many valuable discussions and advises. The work of T. C. is supported 
by the Deutsche Forschungsgemeinschaft in the form of a Forschungsstipendium CH398/1. 
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